Purpose of the script
This script conducts any necessary preprocessing steps to make the data fit for geonization. All result datasets will be in raster format and can be exported as .tif files in the end.
Data
The following data is processed in this scipt.Preparation
Load master mask
Loads a raster that delineates the area of interest. It has a resolution of 1km and is in the right crs. This master mask will be used to format all other datasets to 1km resolution and to clip them to the boudaries of our aoi.
# Loads mask 1km/cell to clip larger datasets to. The mask has exactly th extent and resolution that we aim at.
Mask<-rast("R:\\02_PROJECTS\\01_P_330001\\127_UNCHAIN\\03__Work\\WP4\\03_Casestudy\\02_Data\\Mask_Salzburg_1km\\Mask_Salzburg_1km.tif")
Mask<-aggregate(x=Mask,fact=10, fun=mean)
plot(Mask)
Define common CRS
A crs is defined that can be applied to all datasets that come in with a different crs. It corresponds to the crs of the master mask.
Accessibility
“Accessibility was described as the ease with which a good or service can be reached by the resident population. In a first step, this was measured considering the distance (shortest travel time) along the road network by private car from an origin (resident population) to the nearest facility (infrastructure). This step was implemented in a Geographic Information System using georeferenced input datasets (demographic and infrastructure data) as well as a road network for the calculation of distances. In a second step, the distances to a range of infrastructure facilities were aggregated to indicators of accessibility in an environment for statistical computing.” (R:\02_PROJECTS\01_P_330001\79_RESPECT\03_Work\WP2\Data\Accessibility_STATAT\project_report_accessibility.pdf)
-
CRS [+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs]
-
Cropped
-
Rasterized
-
1km resolution
-
Stretched to 8-bit
- Make N/A 255
- Plot
# Load data
Access<-vect("R:\\02_PROJECTS\\01_P_330001\\79_RESPECT\\03_Work\\WP2\\Data\\Accessibility_STATAT\\Accessibility_Grid.shp")
# Crop Access shp to Mask extent and resolution
Access_rast<-terra::rasterize(Access,Mask,"results__1")
# Stretch to 8-bit
Access_stretch<-stretch(Access_rast, minv=0, maxv=254, smin=0, smax=100)
# Assign N/A to 255
Access_stretch[is.na(Access_stretch)] <- 255
# Visualization
Access_stretch_rast <- raster::raster(Access_stretch)
access <- gplot(Access_stretch_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#b7410e', mid = '#f0e130', high = '#4cbb17', midpoint=128)+
ggtitle("Accessibility")
ggplotly(access, tooltip= "value")Soil function evaluation: Habitat for soil organsism & Natural soil fertility
-
CRS
-
Reclassified [5a → 5, 5b → 6]
-
Cropped
-
Rasterized
-
1km resolution
-
Stretched to 8-bit
-
Take care of N/A
- Plot
Bodenfunktion<-vect("R:\\02_PROJECTS\\01_P_330001\\127_UNCHAIN\\03__Work\\WP4\\03_Casestudy\\02_Data\\Drought_data_Salzburg\\Bodenfunktionsbewertung_Shapefile\\Bodenfunktionsbewertung\\Bodenfunktionsbewertung.shp")
# Reproject data so that it matches the rest of the data
Bodenfkt_reproj <- terra::project(Bodenfunktion,y=crs_all)
# Soil fertility: Reclassify values 5a & 5b to numbers
unique(Bodenfkt_reproj$BTF13b)## [1] "4" "2" "1" "3" "5a" "5b"Bodenfkt_reproj$BTF13b[Bodenfkt_reproj$BTF13b=="5a"]<-5
Bodenfkt_reproj$BTF13b[Bodenfkt_reproj$BTF13b=="5b"]<-6
# Cast characters to integers
Bodenfkt_reproj$BTF13b <- as.integer(Bodenfkt_reproj$BTF13b)
Bodenfkt_reproj$BTF12b <- as.integer(Bodenfkt_reproj$BTF12b)
# Soil organisms----------------
# Rasterize and plot soil organisms
Organisms_rast<-terra::rasterize(Bodenfkt_reproj,Mask,"BTF12b",mean)
Organisms_rast <- cover(Organisms_rast, Mask, values = NA)
Organisms_stretch<-stretch(Organisms_rast, minv=0, maxv=254, smin=3, smax=5)
Organisms_stretch[Organisms_stretch==1]<-0
Organisms_stretch[is.na(Organisms_stretch)] <- 255
# Visualization
Organisms_stretch_rast <- raster(Organisms_stretch)
orga <- gplot(Organisms_stretch_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#b7410e', mid = '#f0e130', high = '#4cbb17', midpoint=128)+
ggtitle("Habitat for soil organisms")
ggplotly(orga, tooltip= "value")# Soil fertility----------------
# Rasterize and plot soil fertility
Fertility_rast<-terra::rasterize(Bodenfkt_reproj,Mask,"BTF13b",mean)
Fertility_rast <- cover ( Fertility_rast, Mask, values = NA)
Fertility_stretch<-stretch(Fertility_rast, minv=0, maxv=254, smin=1, smax=6)
Fertility_stretch[Fertility_stretch==1]<-0
Fertility_stretch[is.na(Fertility_stretch)] <- 255
# Visualization
Fertility_stretch_rast <- raster(Fertility_stretch)
fert <- gplot(Fertility_stretch_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal() + scale_fill_gradient2(low = '#b7410e', mid = '#f0e130', high = '#4cbb17', midpoint=128)+
ggtitle("Natural soil fertility")
ggplotly(fert, tooltip= "value")Soil map: Water conditions & Grassland quality
-
CRS
-
Cropped
-
Rasterized
-
1km resolution
-
Stretched to 8-bit
-
Take care of N/A
- Plot
# Load Bodenkarte as vector
BoKa<-vect("R:\\02_PROJECTS\\01_P_330001\\79_RESPECT\\03_Work\\WP2\\Data\\boka_1km_etrs89_2016\\boka_1km_etrs89_2016\\boka_1km_2016.shp")
# Make raster from vector file, using the extent and resolution from the mask -- That worked really well!
Wasserverhaeltnisse<-terra::rasterize(BoKa,Mask,"wasser")
Wasserverhaeltnisse <- cover(Wasserverhaeltnisse, Mask, values = NA)
Gruenlandwert<-terra::rasterize(BoKa,Mask,"gruenl")
Gruenlandwert <- cover(Gruenlandwert, Mask, values = NA)
# Stretch to 8-bit interval (values range from 0-254; 255 is NA). Smin is 1 because in this dataset 0 holds the value "Nicht beschrieben" which we are not interested in.
Wasserverh_stretch<-stretch(Wasserverhaeltnisse, minv=0, maxv=254, smin=1, smax=9)
Wasserverh_stretch[Wasserverh_stretch==1]<-0
Wasserverh_stretch[is.na(Wasserverh_stretch)] <- 255
Gruenland_stretch<-stretch(Gruenlandwert, minv=0, maxv=254, smin=1, smax=9)
Gruenland_stretch[Gruenland_stretch==1]<-0
Gruenland_stretch[is.na(Gruenland_stretch)] <- 255
# Visualization
# Wasserverhätnisse------
Wasserverh_stretch_rast <- raster(Wasserverh_stretch)
water <- gplot(Wasserverh_stretch_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#b7410e', mid = '#f0e130', high = '#4cbb17', midpoint=128)+
ggtitle("Water conditions")
ggplotly(water, tooltip= "value")#Grünlandwert-------
Gruenland_stretch_rast <- raster(Gruenland_stretch)
gruen <- gplot(Gruenland_stretch_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#b7410e', mid = '#f0e130', high = '#4cbb17', midpoint=128)+
ggtitle("Grassland quality")
ggplotly(gruen, tooltip= "value")Diversity of Plants in Agriculture
The diversity of plants in agriculture and forestry has been used as a sensitivity indicator in the Regional Futures study (Chiari et al. 2013; Lindenthal et al. 2014; Radinger‐Peer et al. 2015): The more diverse cultivation is, the less vulnerable it is, for example to droughts. For calculating the degree of diversity for each 1 km² grid cell, we made use of the Shannon diversity index (SHDI) (Adamczyk and Tiede 2017). The input data was retrieved from INVEKOS. We included the different classes of cultures cultivated on fields to calculate the mentioned diversity index. (R:\02_PROJECTS\01_P_330001\79_RESPECT\03_Work2!Submission-Lucia_M2-2_description_composite_indicators_20180711_LL.pdf)
Anm (Linda): The mentioned input data from INVEKOS could be INVEKOS Feldstücke Österreich 2016 (https://www.data.gv.at/katalog/dataset/fdba4da7-259e-401e-a0ec-b959e5cdf893).
-
CRS
-
Cropped
-
Rasterized
-
1km resolution
-
Stretched to 8-bit
-
Take care of N/A
- Plot
# This raster it taken over as is from the RESPECT projects. Thus, it does not have to be projected.
# Has to be clipped, though.
Div_Plant <- rast("R:\\02_PROJECTS\\01_P_330001\\127_UNCHAIN\\03__Work\\WP4\\03_Casestudy\\02_Data\\DivPlantsAgri.tif")
Div_Plant_crop <- crop(Div_Plant, Mask)
Div_Plant_crop[is.na(Div_Plant_crop)] <- 255
# Visualization
Div_Plant_crop_rast <- raster(Div_Plant_crop)
plant <- gplot(Div_Plant_crop_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#b7410e', mid = '#f0e130', high = '#4cbb17', midpoint=128)+
ggtitle("Diversity of plants in agriculture")
ggplotly(plant, tooltip= "value")Homesteads
-
CRS [+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs]
-
Cropped
-
Rasterized
-
Count Hofstellen per grid cell
-
1km resolution
-
Stretched to 8-bit
-
Take care of N/A
- Plot
# Load
Hofstellen <- vect("R:\\02_PROJECTS\\01_P_330001\\79_RESPECT\\03_Work\\WP2\\Data\\INVEKOS\\invekos_hofstelle_point\\invekos_hofstelle_point.shp")
#Reproject
Hofstellen_reproj <- terra::project(Hofstellen,y=crs_all)
# Crop
Hofstellen_crop <- crop(Hofstellen_reproj, Mask)
# Rasterize
Hofstellen_rast<-terra::rasterize(x = Hofstellen_crop, y = Mask, fun = length)
# To distinguish between NA within and outside of the study area. All NA vaues within the study area are replaced by the value -1 from the mask dataset.
Hofstellen_rast <- cover (Hofstellen_rast, Mask, values = NA)
# Stretch to 8-bit
Hofstellen_stretch<-stretch(Hofstellen_rast, minv=1, maxv=255, smin=1, smax=21)
Hofstellen_stretch[is.na(Hofstellen_stretch)] <- 255
# The -1 from the mask dataset is replaced by the value 0. The function does not recognize the number -1, and sees it instead as a 1. Thus, we replace the value 1 with the value 0 for a correct result.
Hofstellen_stretch[Hofstellen_stretch==1]<-0
# Visualization
Hofstellen_rast_plot <- raster(Hofstellen_stretch)
hoefe <- gplot(Hofstellen_rast_plot) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#4cbb17', mid = '#f0e130', high ='#b7410e' , midpoint=128)+
ggtitle("Homesteads")
ggplotly(hoefe, tooltip= "value")Residents per industry: Primary sector
-
CRS
-
Cropped
-
Rasterized
-
1km resolution
-
Stretched to 8-bit
- Take care of N/A
- Plot
PrimarySector <- vect("R:\\02_PROJECTS\\01_P_330001\\79_RESPECT\\03_Work\\WP2\\Data\\Raster_STATAT\\Wohnbevoelkerung_nach_Wirtschaftszugehoerigkeit_der_Arbeitsstaette.shp")
PrimarySector_reproj <- project(PrimarySector,
y=crs_all)
PrimarySector_rast<-terra::rasterize(PrimarySector_reproj,Mask,"TertiarySe")
PrimarySector_stretch<-stretch(PrimarySector_rast, minv=0, maxv=254, smin=0, smax=3969)
PrimarySector_stretch[is.na(PrimarySector_stretch)] <- 255
# Visualization
PrimarySector_stretch_rast <- raster(PrimarySector_stretch)
prim <- gplot(PrimarySector_stretch_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#4cbb17', mid = '#f0e130', high ='#b7410e' , midpoint=128)
ggtitle("Residents per industry: Primary")## $title
## [1] "Residents per industry: Primary"
##
## attr(,"class")
## [1] "labels"Accomodation and gastronomy
OENACE_I is not a vliad field name ?!
-
CRS
-
Cropped
-
Rasterized
-
1km resolution
-
Stretched to 8-bit
- Take care of N/A
- Plot
BeherbergungGastro <- vect("R:\\02_PROJECTS\\01_P_330001\\79_RESPECT\\03_Work\\WP2\\Data\\Raster_STATAT\\Arbeitsstaetten_nach_Wirtschaftsabschnitt.shp")
BeherbergungGastro_reproj <- project(BeherbergungGastro,
y=crs_all)
# OENACE_I is not a vliad field name ?!
BeherbergungGastro_rast<-terra::rasterize(BeherbergungGastro_reproj,Mask,"OENACE_I")
BeherbergungGastro_stretch<-stretch(BeherbergungGastro_rast, minv=0, maxv=254, smin=0, smax=376)
BeherbergungGastro_stretch[is.na(BeherbergungGastro_stretch)] <- 255
# Visualization
BeherbergungGastro_stretch_rast <- raster(BeherbergungGastro_stretch)
prim <- gplot(BeherbergungGastro_stretch_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#4cbb17', mid = '#f0e130', high ='#b7410e' , midpoint=128)
ggtitle("Accomodation and gastronomy")## $title
## [1] "Accomodation and gastronomy"
##
## attr(,"class")
## [1] "labels"Imperviousness change 2006-2012
-
CRS [+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs]
-
Cropped
-
Rasterized
-
1km resolution
-
Reclassify
-
Stretched to 8-bit
-
Take care of N/A
- Plot
Imperviousness <- rast("R:\\02_PROJECTS\\01_P_330001\\127_UNCHAIN\\03__Work\\WP4\\03_Casestudy\\02_Data\\Imperviousness\\IMC_0612_020m_eu_03035_d03_full.tif")
Imperviousness_reproj <- project(Imperviousness,Mask)
impervious_crop <- crop(Imperviousness_reproj, Mask)
impervious_crop <- cover (Mask,impervious_crop, values = -1)
#"is","becomes""
#100 = unverändert versiegelt
#201 = unverändert unversiegelt
m <- c(100,255, 201,0)
# Make a 2 column matrix from the reclassified vector (first column = "is", second column = "becomes")
rclmat<-matrix(m,ncol=2,byrow=TRUE)
# Reclassify the original values.
imperviousness_crop_reclass<-classify(impervious_crop, rcl=rclmat, include.lowest=FALSE)
# Make N/A values 255
imperviousness_crop_reclass[is.na(imperviousness_crop_reclass)] <- 255
# Visualization
imperviousness_crop_reclass_rast <- raster(imperviousness_crop_reclass)
imper <- gplot(imperviousness_crop_reclass_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#4cbb17', mid = '#f0e130', high ='#b7410e' , midpoint=128)
ggtitle("Imperviousness change")## $title
## [1] "Imperviousness change"
##
## attr(,"class")
## [1] "labels"Agriculturally disadvantaged areas
-
CRS [+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs]
-
Cropped
-
Rasterized
-
1km resolution
-
Stretched to 8-bit
-
Take care of N/A
- Plot
#This one is a bit tricky, because it is nominal data. It is interesting to have it here, in order to distinguish between different structural areas. However it cannot be ranked like the other datasets.
Ben_Geb <- vect("R:\\02_PROJECTS\\01_P_330001\\127_UNCHAIN\\03__Work\\WP4\\03_Casestudy\\02_Data\\Benachteiligte_Gebiete\\BenGeb_6kat_2019_LAEA.shp")
# Reproject data so that it matches the rest of the data
Ben_Geb_reproj <- project(Ben_Geb,
y=crs_all)
# 3 - Berggebiet / Naturräumlicher Teil im Berggebiet
# 4 - Sonst. benachteiligte Gebiete
# 5 - Kleines Gebiet
# 6 - Übergangsgebiet
Benachteiligte_Gebiete<-terra::rasterize(Ben_Geb_reproj,Mask,"BENGEBCODE")
Benachteiligte_Gebiete <- cover (Benachteiligte_Gebiete, Mask, values = NA)
Benachteiligte_Gebiete[Benachteiligte_Gebiete==1]<-0
Benachteiligte_Gebiete_stretch<-stretch(Benachteiligte_Gebiete, minv=0, maxv=254, smin=0, smax=6)
Benachteiligte_Gebiete_stretch[is.na(Benachteiligte_Gebiete_stretch)] <- 255
# Visualization
Benachteiligte_Gebiete_stretch_rast <- raster(Benachteiligte_Gebiete_stretch)
benach <- gplot(Benachteiligte_Gebiete_stretch_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#4cbb17', mid = '#f0e130', high ='#b7410e' , midpoint=128)
ggtitle("Agriculturally disadvantaged areas ")## $title
## [1] "Agriculturally disadvantaged areas "
##
## attr(,"class")
## [1] "labels"Small Woody Features
-
CRS [+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs]
-
Cropped
-
Rasterized
-
1km resolution
-
Reclassify
-
Stretched to 8-bit
-
Take care of N/A
- Plot
Information: The reprojection gave erroneous results. Also, it took very long to finish. Workaround: In ArcGIS Pro, the tool “Mosaic to new raster” was used to combine the two rasters. The same tool also provided the opportunity to choose a crs, so it does not have to be reprojected here. The mosaiced raster as it came from ArcGIS Pro has a gap in it.
#Load the already mosaiced and projected swf dataset
swf<-rast("R:\\02_PROJECTS\\01_P_330001\\127_UNCHAIN\\03__Work\\WP4\\03_Casestudy\\02_Data\\SmallWoodyFeatures\\swf_mosaic_reproj.tif")
# Crop to mask extent
swf_crop<-crop(swf, Mask)
# The reclassify function is also implemented in the terra package, same parameters
# Make a vector with the reclassification values "from", "to", "becomes". We are only interested in the original values 1 and 3, which are being reclassified to 1. The rest of the values are reclassified to 0.
m <- c(0,0,0, 1,1,1, 2,2,0, 3,3,1, 4,255,0)
# Make a 3 column matrix from the reclassified vector
rclmat<-matrix(m,ncol=3,byrow=TRUE)
# Reclassify the original values.
swf_crop_reclass<-classify(swf_crop, rcl=rclmat, include.lowest=FALSE)
#The Small Woody Features dataset has an original resolution of 5m. Therefore it is here aggregated with the factor 200.
swf_agg1km<-aggregate(x=swf_crop_reclass,fact=200, fun=mean)
# Stretch to 8-bit
swf_stretch<-stretch(swf_agg1km, minv=0, maxv=254, smin=0, smax=1)
swf_stretch[is.na(swf_stretch)] <- 255
# Visualization
swf_agg1km_rast <- raster(swf_stretch)
swf <- gplot(swf_agg1km_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#b7410e', mid = '#f0e130', high ='#4cbb17' , midpoint=128)
ggtitle("Small woody features")## $title
## [1] "Small woody features"
##
## attr(,"class")
## [1] "labels"Landscape- and lake-conservation area
Gebiete außerhalb geschlossener Ortschaften von besonderer landschaftlicher Schönheit und/oder von Bedeutung für die Erholung als charakteristische Naturlandschaft oder naturnahe Kulturlandschaft. Die Ausweisung erfolgt durch Verordnung der Landesregierung. In allen Landschaftsschutzgebieten gilt die Allgemeine Landschaftsschutzverordnung (ALV). Bestimmte Maßnahmen sind nur mit Bewilligung der Bezirksverwaltungsbehörde zulässig. (https://www.salzburg.gv.at/themen/natur/naturschutzrecht-2/naturschutzrecht-salzburg/gebietsschutz/landschaftsschutzgebiet)
-
CRS
-
Cropped
-
Rasterized
-
1km resolution
-
Stretched to 8-bit
- Take care of N/A
- Plot
255 <- not within AOI
254 <- Landschaftsschutzgebiete
0 <- within AOI but not a Landschaftsschutzgebiet
Landschaftsschutz <- vect("R:\\02_PROJECTS\\01_P_330001\\127_UNCHAIN\\03__Work\\WP4\\03_Casestudy\\02_Data\\Naturschutz\\Landschaftsschutzgebiete_Shapefile\\Landschaftsschutzgebiete\\Landschaftsschutzgebiete.shp")
Landschaftsschutz_reproj <- project(Landschaftsschutz,
y=crs_all)
Landschaftsschutz_reproj$yes <- 20
# OENACE_I is not a vliad field name ?
Landschaftsschutz_rast<-terra::rasterize(Landschaftsschutz_reproj,Mask,"yes")
Landschaftsschutz_rast <- cover(Landschaftsschutz_rast, Mask, values = NA)
Landschaftsschutz_stretch<-stretch(Landschaftsschutz_rast, minv=0, maxv=254, smin=-1, smax=20)
Landschaftsschutz_stretch[Landschaftsschutz_stretch==1]<-0
Landschaftsschutz_stretch[is.na(Landschaftsschutz_stretch)] <- 255
# Visualization
Landschaftsschutz_stretch_rast <- raster(Landschaftsschutz_stretch)
prim <- gplot(Landschaftsschutz_stretch_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#4cbb17', mid = '#f0e130', high ='#b7410e' , midpoint=128)
ggtitle("Accomodation and gastronomy")## $title
## [1] "Accomodation and gastronomy"
##
## attr(,"class")
## [1] "labels"Nature reserve
Gebiete mit völliger bis weitgehender Ursprünglichkeit und/oder Gebiete, die seltene oder gefährdete Tier- oder Pflanzenarten oder solche Lebensgemeinschaften von Tieren oder Pflanzen aufweisen. Die Ausweisung erfolgt durch Verordnung der Landesregierung. Bestimmte Eingriffe in Naturschutzgebiete sind grundsätzlich verboten oder nur mit Bewilligung der Landesregierung zulässig.(https://www.salzburg.gv.at/themen/natur/naturschutzrecht-2/naturschutzrecht-salzburg/gebietsschutz/naturschutzgebiet)
255 <- not within AOI
254 <- Naturschutzgebiete
0 <- within AOI but not a Naturschutzgebiet
Naturschutz <- vect("R:\\02_PROJECTS\\01_P_330001\\127_UNCHAIN\\03__Work\\WP4\\03_Casestudy\\02_Data\\Naturschutz\\Naturschutzgebiete_Shapefile\\Naturschutzgebiete\\Naturschutzgebiete.shp")
Naturschutz_reproj <- project(Naturschutz,
y=crs_all)
Naturschutz_reproj$yes <- 20
# OENACE_I is not a vliad field name ?
Naturschutz_rast<-terra::rasterize(Naturschutz_reproj,Mask,"yes")
Naturschutz_rast <- cover(Naturschutz_rast, Mask, values = NA)
Naturschutz_stretch<-stretch(Naturschutz_rast, minv=0, maxv=254, smin=-1, smax=20)
Naturschutz_stretch<-stretch(Naturschutz_rast, minv=0, maxv=254, smin=0, smax=3)
Naturschutz_stretch[is.na(Naturschutz_stretch)] <- 255
# Visualization
Naturschutz_stretch_rast <- raster(Naturschutz_stretch)
prim <- gplot(Naturschutz_stretch_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#4cbb17', mid = '#f0e130', high ='#b7410e' , midpoint=128)
ggtitle("Nature reserve")## $title
## [1] "Nature reserve"
##
## attr(,"class")
## [1] "labels"Forested area
255 <- Forest outside AOI
254 <- Forested area
0 <- Non forested area
Forest <- vect("R:\\02_PROJECTS\\01_P_330001\\127_UNCHAIN\\03__Work\\WP4\\03_Casestudy\\02_Data\\Drought_data_Salzburg\\Waldflaechen_Shapefile\\Waldflaechen\\Waldflaechen.shp")
Forest_reproj <- project(Forest,
y=crs_all)
Forest_reproj$yes <- 20
Forest_rast<-terra::rasterize(Forest_reproj,Mask,"yes")
Forest_rast <- cover(Forest_rast, Mask, values = NA)
Forest_stretch<-stretch(Forest_rast, minv=0, maxv=254, smin=-1, smax=20)
Forest_stretch[is.na(Forest_stretch)] <- 255
# Visualization
Forest_stretch_rast <- raster(Forest_stretch)
prim <- gplot(Forest_stretch_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#4cbb17', mid = '#f0e130', high ='#b7410e' , midpoint=128)
ggtitle("Forested area")## $title
## [1] "Forested area"
##
## attr(,"class")
## [1] "labels"Industrial water supply
Nutzwasserversorgung & Trink-und Nutzwasserversorgung für betriebliche Versorgung.
Wasser <- vect("R:\\02_PROJECTS\\01_P_330001\\127_UNCHAIN\\03__Work\\WP4\\03_Casestudy\\02_Data\\alle_Wasserpunkte\\alle_Wasserpunkte.shp")
Wasser_nutzwasser <- subset(Wasser, Wasser$SPRT_NAME %in% c ("Nutzwasserversorgung", "Trink-und Nutzwasserversorgung für betriebliche Versorgung"))
Wasser_reproj <- project(Wasser_nutzwasser,
y=crs_all)
Wasser_reproj$yes <- 20
Wasser_rast<-terra::rasterize(Wasser_reproj, Mask, fun = length)
Wasser_rast <- cover(Wasser_rast, Mask, values = NA)
Wasser_stretch<-stretch(Wasser_rast, minv=0, maxv=254, smin=-1, smax=20)
Wasser_stretch[is.na(Wasser_stretch)] <- 255
Wasser_stretch[Wasser_stretch==1]<-0
# Visualization
Wasser_stretch_rast <- raster(Wasser_stretch)
prim <- gplot(Wasser_stretch_rast) + geom_tile(aes(fill=value), alpha=0.8) + coord_equal()+ scale_fill_gradient2(low = '#4cbb17', mid = '#f0e130', high ='#b7410e' , midpoint=128)
ggtitle("Nutzwasserversorgung")## $title
## [1] "Nutzwasserversorgung"
##
## attr(,"class")
## [1] "labels"Export
[x] Accessibility
[x] Soil function evaluation: Habitat for soil organsism
[x] Soil function evaluation: Natural soil fertility
[x] Soil map: Water conditions
[x] Soil map: Gassland quality
[x] Diversity of Plants in Agriculture
[x] Homesteads
[x] Imperviousness change
[x] Agriculturally diadvantaged areas
[x] Verdichtungsgefährdete Flächen
[x] Small Woody Features
x Standardized Precipitation Index (SPI)
[x] Inhabitants per sector: Primary sector
final_data_path <- "R:\\02_PROJECTS\\01_P_330001\\127_UNCHAIN\\03__Work\\WP4\\03_Casestudy\\02_Data\\00_Final_data"
data_cube <- c(Access_stretch, swf_stretch, Organisms_stretch, Fertility_stretch, Wasserverh_stretch, Gruenland_stretch, Div_Plant_crop, Hofstellen_rast, PrimarySector_stretch, imperviousness_crop_reclass, Benachteiligte_Gebiete_stretch)
# Accessibility
ac <- file.path(final_data_path,"accessibility.tif")
terra::writeRaster(Access_stretch, ac, overwrite = TRUE)
ac <- file.path(final_data_path,"accessibility.tif")
terra::writeRaster(Access_stretch, ac, overwrite = TRUE)
f <- file.path(final_data_path,"small_woody_features.tif")
terra::writeRaster(swf_stretch, f, overwrite = TRUE)
f <- file.path(final_data_path,"data_cube.tif")
terra::writeRaster(data_cube, f, overwrite = TRUE)
#library(stars)
#library(gdalUtils)
#imperviousness_crop_reclass %>%
# st_as_stars %>%
# write_stars("Imperviousness.gpkg",
# driver = "GPKG")
#swf_crop %>%
# st_as_stars %>%
# write_stars("SmallWoodyFeatures.gpkg",
# driver = "GPKG")
#swf_agg1km %>%
# st_as_stars %>%
# write_stars("SmallWoodyFeatures1km.gpkg",
# driver = "GPKG")
#plough_crop %>%
# st_as_stars %>%
# write_stars("Plughing_6years.gpkg",
# driver = "GPKG")
#gdalUtils::gdalinfo("Drought_casestudy_rasters.gpkg") %>%
# cat(sep = "\n")